rm(list = ls())
set.seed(12345678)
library(bpCausal)
library(modeest)



## parallel computing
#library(doParallel)
#library(doRNG)
#library(future)


## setwd("/Users/lichengl/Desktop/BSA_replication_code/")

source("code/sim.R")
source("code/plot_function.R")
library(parallel)

## 500 simulations 
nboots <- 1000





## ---------------------- ##
## case 1: N = 50, T = 30 ##
## ---------------------- ##


## save the Monte Carlo results: 4 is number of models, 10000 is mcmc iternation numbers 
mc <- array(NA, dim = c(4, 10000, nboots))


N <- 50   ## 50, 100
TT <- 30  ## 30, 50, 90 
Ntr <- N * 0.2 ## 10, 20
T0 <- TT - 10  ## 20, 40, 80 
p <- 5
beta <- c(1, 3, 0, 0, 0)
## beta <- c(1, 3)
r <- 2
twoway <- 1
c12 <- 0.25
c22 <- 0.25
c32 <- 0.25
c42 <- 0.25
c52 <- 0.25


Xname <- c("X1", "X2", "X3", "X4", "X5")
Xname2 <- c("X1", "X2", "U") ## oracal case where we know U 


## with one unobserevd confounder 
## no unobsevred confounder 

oneboot <- function(kk) {

	datas <- simulate(N = N, TT = TT, Ntr = Ntr, T0 = T0, p = p,
		           beta = beta, r = 2, twoway = twoway, ATT = 0, 
		           c12 = c12, c22 = c22, c32 = c32,
		           c42 = c42, c52 = c52, coeff = 0.6, betau = 1) ## 1 unobserved confounder, sharp null case, it doesn't matter...

	data <- datas$data 

	## mean(data$U[which(data$D == 1)]) - mean(data$U[which(data$D == 0)]) 

	sub_ff <- matrix(NA, 4, 10000)



	## proposed model
	fit1 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 1, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)



	## sense with diffuse prior
	fit2 <-  bpCausal(data = data, index = c("id", "time"), 
		              Yname = "Y", Dname = "D", Xname = Xname, 
		              Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25, sigmau2 = 10)


	## naive case 
	fit0 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## ideal case
	fit3 <-  bpCausal(data = data, index = c("id", "time"), 
	 	           Yname = "Y", Dname = "D", Xname = Xname2, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## mean(datas$veff[(T0+1):TT])

	avg2 <- effSummary_simple(fit2)$est.avg
	avg1 <- effSummary_simple(fit1)$est.avg
	avg0 <- effSummary_simple(fit0)$est.avg
	avg3 <- effSummary_simple(fit3)$est.avg



    sub_ff[1, ] <- avg0 
    sub_ff[2, ] <- avg1 
    sub_ff[3, ] <- avg2 
    sub_ff[4, ] <- avg3 

	rm(fit2, fit1, fit0, fit3, 
	   avg2, avg1, avg0, avg3)
	
	gc(TRUE) ## release memory

	
    return(sub_ff)

}


## I use 10 core parallels 

for (i in 1:100) {	
    
    pos <- (i - 1) * 10
	subf <- mclapply(1:10 + pos, oneboot, mc.cores = 10)

	for (j in 1:10) {
		
		mc[,,pos + j] <- subf[[j]]
	}
	cat(paste("iterations: ", i, sep = ""))
	cat("\n")
}


## rename and save results
mc_50_30 <- mc 

save(mc_50_30, file = "results/sim_N50_T30.RData")




## ---------------------- ##
## case 2: N = 50, T = 50 ##
## ---------------------- ##


## save the Monte Carlo results: 4 is number of models, 10000 is mcmc iternation numbers 
mc <- array(NA, dim = c(4, 10000, nboots))


N <- 50   ## 50, 100
TT <- 50  ## 30, 50, 90 
Ntr <- N * 0.2 ## 10, 20
T0 <- TT - 10  ## 20, 40, 80 
p <- 5
beta <- c(1, 3, 0, 0, 0)
## beta <- c(1, 3)
r <- 2
twoway <- 1
c12 <- 0.25
c22 <- 0.25
c32 <- 0.25
c42 <- 0.25
c52 <- 0.25


Xname <- c("X1", "X2", "X3", "X4", "X5")
Xname2 <- c("X1", "X2", "U") ## oracal case where we know U 


## with one unobserevd confounder 
## no unobsevred confounder 

oneboot <- function(kk) {

	datas <- simulate(N = N, TT = TT, Ntr = Ntr, T0 = T0, p = p,
		           beta = beta, r = 2, twoway = twoway, ATT = 0, 
		           c12 = c12, c22 = c22, c32 = c32,
		           c42 = c42, c52 = c52, coeff = 0.6, betau = 1) ## 1 unobserved confounder, sharp null case, it doesn't matter...

	data <- datas$data 

	## mean(data$U[which(data$D == 1)]) - mean(data$U[which(data$D == 0)]) 

	sub_ff <- matrix(NA, 4, 10000)



	## proposed model
	fit1 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 1, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)



	## sense with diffuse prior
	fit2 <-  bpCausal(data = data, index = c("id", "time"), 
		              Yname = "Y", Dname = "D", Xname = Xname, 
		              Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25, sigmau2 = 10)


	## naive case 
	fit0 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## ideal case
	fit3 <-  bpCausal(data = data, index = c("id", "time"), 
	 	           Yname = "Y", Dname = "D", Xname = Xname2, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## mean(datas$veff[(T0+1):TT])

	avg2 <- effSummary_simple(fit2)$est.avg
	avg1 <- effSummary_simple(fit1)$est.avg
	avg0 <- effSummary_simple(fit0)$est.avg
	avg3 <- effSummary_simple(fit3)$est.avg



    sub_ff[1, ] <- avg0 
    sub_ff[2, ] <- avg1 
    sub_ff[3, ] <- avg2 
    sub_ff[4, ] <- avg3 

	rm(fit2, fit1, fit0, fit3, 
	   avg2, avg1, avg0, avg3)
	
	gc(TRUE) ## release memory

	
    return(sub_ff)

}


## I use 10 core parallels 

for (i in 1:100) {	
    
    pos <- (i - 1) * 10
	subf <- mclapply(1:10 + pos, oneboot, mc.cores = 10)

	for (j in 1:10) {
		
		mc[,,pos + j] <- subf[[j]]
	}
	cat(paste("iterations: ", i, sep = ""))
	cat("\n")
}


## rename and save results
mc_50_50 <- mc 

save(mc_50_50, file = "results/sim_N50_T50.RData")


## ---------------------- ##
## case 3: N = 50, T = 90 ##
## ---------------------- ##


## save the Monte Carlo results: 4 is number of models, 10000 is mcmc iternation numbers 
mc <- array(NA, dim = c(4, 10000, nboots))


N <- 50   ## 50, 100
TT <- 90  ## 30, 50, 90 
Ntr <- N * 0.2 ## 10, 20
T0 <- TT - 10  ## 20, 40, 80 
p <- 5
beta <- c(1, 3, 0, 0, 0)
## beta <- c(1, 3)
r <- 2
twoway <- 1
c12 <- 0.25
c22 <- 0.25
c32 <- 0.25
c42 <- 0.25
c52 <- 0.25


Xname <- c("X1", "X2", "X3", "X4", "X5")
Xname2 <- c("X1", "X2", "U") ## oracal case where we know U 


## with one unobserevd confounder 
## no unobsevred confounder 

oneboot <- function(kk) {

	datas <- simulate(N = N, TT = TT, Ntr = Ntr, T0 = T0, p = p,
		           beta = beta, r = 2, twoway = twoway, ATT = 0, 
		           c12 = c12, c22 = c22, c32 = c32,
		           c42 = c42, c52 = c52, coeff = 0.6, betau = 1) ## 1 unobserved confounder, sharp null case, it doesn't matter...

	data <- datas$data 

	## mean(data$U[which(data$D == 1)]) - mean(data$U[which(data$D == 0)]) 

	sub_ff <- matrix(NA, 4, 10000)



	## proposed model
	fit1 <- bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 1, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)



	## sense with diffuse prior
	fit2 <- bpCausal(data = data, index = c("id", "time"), 
		              Yname = "Y", Dname = "D", Xname = Xname, 
		              Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25, sigmau2 = 10)


	## naive case 
	fit0 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## ideal case
	fit3 <-  bpCausal(data = data, index = c("id", "time"), 
	 	           Yname = "Y", Dname = "D", Xname = Xname2, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## mean(datas$veff[(T0+1):TT])

	avg2 <- effSummary_simple(fit2)$est.avg
	avg1 <- effSummary_simple(fit1)$est.avg
	avg0 <- effSummary_simple(fit0)$est.avg
	avg3 <- effSummary_simple(fit3)$est.avg



    sub_ff[1, ] <- avg0 
    sub_ff[2, ] <- avg1 
    sub_ff[3, ] <- avg2 
    sub_ff[4, ] <- avg3 

	rm(fit2, fit1, fit0, fit3, 
	   avg2, avg1, avg0, avg3)
	
	gc(TRUE) ## release memory

	
    return(sub_ff)

}


## I use 10 core parallels 

for (i in 1:100) {	
    
    pos <- (i - 1) * 10
	subf <- mclapply(1:10 + pos, oneboot, mc.cores = 10)

	for (j in 1:10) {
		
		mc[,,pos + j] <- subf[[j]]
	}
	cat(paste("iterations: ", i, sep = ""))
	cat("\n")
}


## rename and save results
mc_50_90 <- mc 

save(mc_50_90, file = "results/sim_N50_T90.RData")




## ----------------------- ##
## case 4: N = 100, T = 30 ##
## ----------------------- ##


## save the Monte Carlo results: 4 is number of models, 10000 is mcmc iternation numbers 
mc <- array(NA, dim = c(4, 10000, nboots))


N <- 100   ## 50, 100
TT <- 30  ## 30, 50, 90 
Ntr <- N * 0.2 ## 10, 20
T0 <- TT - 10  ## 20, 40, 80 
p <- 5
beta <- c(1, 3, 0, 0, 0)
## beta <- c(1, 3)
r <- 2
twoway <- 1
c12 <- 0.25
c22 <- 0.25
c32 <- 0.25
c42 <- 0.25
c52 <- 0.25


Xname <- c("X1", "X2", "X3", "X4", "X5")
Xname2 <- c("X1", "X2", "U") ## oracal case where we know U 


## with one unobserevd confounder 
## no unobsevred confounder 

oneboot <- function(kk) {

	datas <- simulate(N = N, TT = TT, Ntr = Ntr, T0 = T0, p = p,
		           beta = beta, r = 2, twoway = twoway, ATT = 0, 
		           c12 = c12, c22 = c22, c32 = c32,
		           c42 = c42, c52 = c52, coeff = 0.6, betau = 1) ## 1 unobserved confounder, sharp null case, it doesn't matter...

	data <- datas$data 

	## mean(data$U[which(data$D == 1)]) - mean(data$U[which(data$D == 0)]) 

	sub_ff <- matrix(NA, 4, 10000)



	## proposed model
	fit1 <- bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 1, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)



	## sense with diffuse prior
	fit2 <- bpCausal(data = data, index = c("id", "time"), 
		              Yname = "Y", Dname = "D", Xname = Xname, 
		              Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25, sigmau2 = 10)


	## naive case 
	fit0 <- bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## ideal case
	fit3 <- bpCausal(data = data, index = c("id", "time"), 
	 	           Yname = "Y", Dname = "D", Xname = Xname2, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## mean(datas$veff[(T0+1):TT])

	avg2 <- effSummary_simple(fit2)$est.avg
	avg1 <- effSummary_simple(fit1)$est.avg
	avg0 <- effSummary_simple(fit0)$est.avg
	avg3 <- effSummary_simple(fit3)$est.avg



    sub_ff[1, ] <- avg0 
    sub_ff[2, ] <- avg1 
    sub_ff[3, ] <- avg2 
    sub_ff[4, ] <- avg3 

	rm(fit2, fit1, fit0, fit3, 
	   avg2, avg1, avg0, avg3)
	
	gc(TRUE) ## release memory

	
    return(sub_ff)

}


## I use 10 core parallels 

for (i in 1:100) {	
    
    pos <- (i - 1) * 10
	subf <- mclapply(1:10 + pos, oneboot, mc.cores = 10)

	for (j in 1:10) {
		
		mc[,,pos + j] <- subf[[j]]
	}
	cat(paste("iterations: ", i, sep = ""))
	cat("\n")
}


## rename and save results
mc_100_30 <- mc 

save(mc_100_30, file = "results/sim_N100_T30.RData")





## ----------------------- ##
## case 5: N = 100, T = 50 ##
## ----------------------- ##


## save the Monte Carlo results: 4 is number of models, 10000 is mcmc iternation numbers 
mc <- array(NA, dim = c(4, 10000, nboots))


N <- 100   ## 50, 100
TT <- 50  ## 30, 50, 90 
Ntr <- N * 0.2 ## 10, 20
T0 <- TT - 10  ## 20, 40, 80 
p <- 5
beta <- c(1, 3, 0, 0, 0)
## beta <- c(1, 3)
r <- 2
twoway <- 1
c12 <- 0.25
c22 <- 0.25
c32 <- 0.25
c42 <- 0.25
c52 <- 0.25


Xname <- c("X1", "X2", "X3", "X4", "X5")
Xname2 <- c("X1", "X2", "U") ## oracal case where we know U 


## with one unobserevd confounder 
## no unobsevred confounder 

oneboot <- function(kk) {

	datas <- simulate(N = N, TT = TT, Ntr = Ntr, T0 = T0, p = p,
		           beta = beta, r = 2, twoway = twoway, ATT = 0, 
		           c12 = c12, c22 = c22, c32 = c32,
		           c42 = c42, c52 = c52, coeff = 0.6, betau = 1) ## 1 unobserved confounder, sharp null case, it doesn't matter...

	data <- datas$data 

	## mean(data$U[which(data$D == 1)]) - mean(data$U[which(data$D == 0)]) 

	sub_ff <- matrix(NA, 4, 10000)



	## proposed model
	fit1 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 1, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)



	## sense with diffuse prior
	fit2 <-  bpCausal(data = data, index = c("id", "time"), 
		              Yname = "Y", Dname = "D", Xname = Xname, 
		              Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25, sigmau2 = 10)


	## naive case 
	fit0 <-  bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## ideal case
	fit3 <-  bpCausal(data = data, index = c("id", "time"), 
	 	           Yname = "Y", Dname = "D", Xname = Xname2, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## mean(datas$veff[(T0+1):TT])

	avg2 <- effSummary_simple(fit2)$est.avg
	avg1 <- effSummary_simple(fit1)$est.avg
	avg0 <- effSummary_simple(fit0)$est.avg
	avg3 <- effSummary_simple(fit3)$est.avg



    sub_ff[1, ] <- avg0 
    sub_ff[2, ] <- avg1 
    sub_ff[3, ] <- avg2 
    sub_ff[4, ] <- avg3 

	rm(fit2, fit1, fit0, fit3, 
	   avg2, avg1, avg0, avg3)
	
	gc(TRUE) ## release memory

	
    return(sub_ff)

}


## I use 10 core parallels 

for (i in 1:100) {	
    
    pos <- (i - 1) * 10
	subf <- mclapply(1:10 + pos, oneboot, mc.cores = 10)

	for (j in 1:10) {
		
		mc[,,pos + j] <- subf[[j]]
	}
	cat(paste("iterations: ", i, sep = ""))
	cat("\n")
}


## rename and save results
mc_100_50 <- mc 

save(mc_100_50, file = "results/sim_N100_T50.RData")





## ----------------------- ##
## case 6: N = 100, T = 90 ##
## ----------------------- ##


## save the Monte Carlo results: 4 is number of models, 10000 is mcmc iternation numbers 
mc <- array(NA, dim = c(4, 10000, nboots))


N <- 100   ## 50, 100
TT <- 90  ## 30, 50, 90 
Ntr <- N * 0.2 ## 10, 20
T0 <- TT - 10  ## 20, 40, 80 
p <- 5
beta <- c(1, 3, 0, 0, 0)
## beta <- c(1, 3)
r <- 2
twoway <- 1
c12 <- 0.25
c22 <- 0.25
c32 <- 0.25
c42 <- 0.25
c52 <- 0.25


Xname <- c("X1", "X2", "X3", "X4", "X5")
Xname2 <- c("X1", "X2", "U") ## oracal case where we know U 


## with one unobserevd confounder 
## no unobsevred confounder 

oneboot <- function(kk) {

	datas <- simulate(N = N, TT = TT, Ntr = Ntr, T0 = T0, p = p,
		           beta = beta, r = 2, twoway = twoway, ATT = 0, 
		           c12 = c12, c22 = c22, c32 = c32,
		           c42 = c42, c52 = c52, coeff = 0.6, betau = 1) ## 1 unobserved confounder, sharp null case, it doesn't matter...

	data <- datas$data 

	## mean(data$U[which(data$D == 1)]) - mean(data$U[which(data$D == 0)]) 

	sub_ff <- matrix(NA, 4, 10000)



	## proposed model
	fit1 <- bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 1, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)



	## sense with diffuse prior
	fit2 <- bpCausal(data = data, index = c("id", "time"), 
		              Yname = "Y", Dname = "D", Xname = Xname, 
		              Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 1.001, p2 = 0.001, 
	                  sense = 1,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25, sigmau2 = 10)


	## naive case 
	fit0 <- bpCausal(data = data, index = c("id", "time"), 
		           Yname = "Y", Dname = "D", Xname = Xname, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## ideal case
	fit3 <- bpCausal(data = data, index = c("id", "time"), 
	 	           Yname = "Y", Dname = "D", Xname = Xname2, 
		           Zname = NULL, Aname = NULL, re = "both", 
	                  ar1 = FALSE, r = 2, niter = 11000, burn = 1000, 
	                  xlasso = 0, zlasso = 0, alasso = 0, flasso = 0, 
	                  a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001, 
	                  c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001, 
	                  sense = 0,  c12 = 0.25, c22 = 0.25, 
	                  c32 = 0.25, c42 = 0.25, c52 = 0.25)


	## mean(datas$veff[(T0+1):TT])

	avg2 <- effSummary_simple(fit2)$est.avg
	avg1 <- effSummary_simple(fit1)$est.avg
	avg0 <- effSummary_simple(fit0)$est.avg
	avg3 <- effSummary_simple(fit3)$est.avg



    sub_ff[1, ] <- avg0 
    sub_ff[2, ] <- avg1 
    sub_ff[3, ] <- avg2 
    sub_ff[4, ] <- avg3 

	rm(fit2, fit1, fit0, fit3, 
	   avg2, avg1, avg0, avg3)
	
	gc(TRUE) ## release memory

	
    return(sub_ff)

}


## I use 10 core parallels 

for (i in 1:100) {	
    
    pos <- (i - 1) * 10
	subf <- mclapply(1:10 + pos, oneboot, mc.cores = 10)

	for (j in 1:10) {
		
		mc[,,pos + j] <- subf[[j]]
	}
	cat(paste("iterations: ", i, sep = ""))
	cat("\n")
}


## rename and save results
mc_100_90 <- mc 

save(mc_100_90, file = "results/sim_N100_T90.RData")


